function[u,v,xb,yb,ub,vb]=particle_coord(n,R,t,om,om1,om2,phi,mu,epsilon,xout,yout,Nout,dt0)

% COORDINATES AND VELOCITY OF THE MOVING BOUNDARY
phi1=0;
modom_old=0;
modom=mod(floor((om*t+phi)/pi),2);
if modom==0
    if modom_old==1
        % ANY TIME modom SWITCHES FROM 1 TO 0 THE PHASE IS REDEFINED
        phi1=om2*(t-dt0)+phi2-om1*(t-dt0);
    end
        [xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,om1,phi1);
        om=om1;
        phi=phi1;
else
    if modom_old==0
        % ANY TIME modom SWITCHES FROM 0 TO 1 THE PHASE IS REDEFINED
        phi2=om1*(t-dt0)+phi1-om2*(t-dt0);
    end
    [xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,om2,phi2);
    om=om2;
    phi=phi2;
end
modom_old=modom;
% clf;
% figure(99)
% plot(xb, yb, 'k.')
% title(['t= ' num2str(t)])
% xlabel('X')
% ylabel('Y')
% daspect([1 1 1])
% axis([-2*R 2*R -2*R 2*R])
% grid on
% pause(.1)
% A(:,i) = getframe(gcf);
% movie2avi(A, 'boundary.avi', 'fps', 5, 'quality', 100);
% THE REGULARIZED GREEN FUNCTION MATRIX IS COMPUTED
[G]=G_matrix_mod2D(xb,yb,n,xb,yb,n,mu,epsilon);
% figure
% imagesc(G)
% colorbar
% title('G')
for j=1:n
    vel_in((j-1)*2+1)=ub(j);
	vel_in((j-1)*2+2)=vb(j);
end
% SOLVE THE LINEAR SYSTEM TO FIND THE FORCING f ON THE BOUNDARY
opts.SYM = true;
f=linsolve(G,vel_in',opts);
[M2N]=G_matrix_mod2D(xout,yout,Nout,xb,yb,n,mu,epsilon);
uout=M2N*f;
for j=1:Nout
    u(j)=uout((j-1)*2+1);
	v(j)=uout((j-1)*2+2);
end